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THE PIECEWISE PARABOLIC METHOD FOR 
MULTIDIMENSIONAL RELATIVISTIC FLUID DYNAMICS 

A. Mignone 1 - 2 ' 3 , T. Plewa 2 ' 3 ' 4 , and G. Bode- 1 
^ ■ ABSTRACT 

o 

' We present an extension of the Piecewise Parabolic Method to special relativistic fluid dy- 

namics in multidimcnsions. The scheme is conservative, dimensionally unsplit, and suitable for a 
general equation of state. Temporal evolution is second-order accurate and employs characteristic 
projection operators; spatial interpolation is piece-wise parabolic making the scheme third-order 
accurate in smooth regions of the flow away from discontinuities. The algorithm is written for 
a general system of orthogonal curvilinear coordinates and can be used for computations in 
non-cartesian geometries. A non-linear iterative Riemann solver based on the two-shock approx- 
■ imation is used in flux calculation. In this approximation, an initial discontinuity decays into 

a set of discontinuous waves only implying that, in particular, rarefaction waves are treated as 
flow discontinuities. We also present a new and simple equation of state which approximates the 
ly-j . exact result for the relativistic perfect gas with high accuracy. The strength of the new method is 

demonstrated in a series of numerical tests and more complex simulations in one, two and three 
*~Q ■ dimensions. 

o 

' Subject headings: hydrodynamics — relativity — shock waves — methods: numerical 
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INTRODUCTION 



Highly energetic astrophysical phenomena are known to be, in many cases, relativistic in nature. A wide 
range of objects, in fact, exhibits a number of properties that can be accounted for only in the framework 
of the theory of special or general relativity: superluminal motion of relativistic jets in extragalactic radio 
sources (Begelman et al. 1984), jets and accretion flows around massive compact objects (Koide et al. 1999; 
Meier et al. 2001), pulsar winds (Del Zanna et al. 2004; Bogovalov et al. 2005), gamma ray bursts (Aloy et 
al. 2000; Zhang et al. 2003; Mizuno et al. 2004), as well as particle beams produced in heavy- ion collisions in 
terrestrial experiments (Ackermann et al. 2001; Morita et al. 2002; Molnar & Huovinen 2004) characterized by 
flow velocities very close (> 0.99c) to the speed of light. Such problems are difficult to study numerically due 
to inherent complexity of the problem itself exhibiting rapid spatial and temporal changes often demanding 
development of more efficient algorithms. 

Recent progress in the development of powerful methods in computational fluid dynamics allowed re- 
searchers to gain understanding of relativistic flows. Numerical formulations of the special relativistic fluid 
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equations have been investigated by several authors in the last decade (for an excellent review see Marti & 
Miiller 2003). However, there now exists strong evidence that a particular class of numerical methods, the 
so-called high-resolution shock-capturing schemes (HRSC henceforth, e.g., Donat et al. 1998; Marquina et 
al. 1992), provide the necessary means to develop stable and robust relativistic (special or general) fluid dy- 
namics codes. HRSC methods rely on the conservative formulation of the fluid equations. This formulation 
is of fundamental importance in representing the evolution of flows with steep gradients and discontinuities. 
One of the key aspects of these schemes is the temporal evolution of the system, which routinely involves 
(exact or approximate) solution of Riemann problems at the interfaces separating numerical grid cells, thus 
making these schemes highly successful in modeling discontinuous flows. Besides, these algorithms have at 
least 2nd order accuracy in smooth regions of the flow and have been shown to produce excellent results for 
a wide variety of problems. 

Shock-capturing schemes have a long standing tradition in the framework of solving the Euler equations 
of gas dynamics (e.g., Colella 1985; Colella & Woodward 1984; VanLeer 1997; Toro 1997). Several of these 
"classical" schemes have now been extended to special relativistic hydrodynamics (see Balsara 1994; Dai & 
Woodward 1997; Falle & Komissarov 1996; Donat et al. 1998; Sokolov et al. 2001, and references therein). The 
Piecewise Parabolic Method (PPM) by Colella & Woodward (1984, CW84 henceforth) is still considered as 
the state-of-the-art method in computational fluid dynamics. The original PPM algorithm has been recently 
re-formulated for the relativistic fluid flows in one spatial dimension by Marti & Miiller (1996). Here we 
present an extension of the PPM method to multidimensional relativistic hydrodynamics. As wc will show, 
the main differences between the classical and relativistic version of PPM are the coupling between normal 
and transverse velocity components introduced by the Lorentz factor and the coupling of the latter to the 
specific enthalpy. 

Our relativistic PPM consists of several components. The piece-wise parabolic interpolation is done 
in the volume coordinate (CW84) and provides third-order accuracy in space in smooth parts of the flow. 
Temporal evolution uses characteristic information and is second-order accurate. The scheme also uses a 
non-linear iterative Riemann solver based on the two-shock approximation. The unsplit fully-coupled corner- 
transport upwind method is used for advection (Colella 1990; Saltzman 1994). Several verifications tests are 
presented to demonstrate the correctness of the implementation. 

The paper is organized as follows. In §2 we give a short review of the special relativistic fluid equations. 
In §3 we present discretization of the equations in a general system of orthogonal curvilinear coordinates and 
describe the numerical method. We include the details of the time evolution algorithm based on characteristic 
tracing. The Riemann solver is presented in §3.3. In §3.4 we consider different equation of states and 
introduce a new formulation suitable for relativistic regimes. Several numerical tests are presented in §4 and 
conclusions are drawn in §5. 



2. THE EQUATIONS OF RELATIVISTIC HYDRODYNAMICS 

In the framework of special relativity, the motion of an ideal fluid is governed by the laws of particle 
number conservation and energy-momentum conservation (Landau & Lifshitz 1959; Weinberg 1972). In the 
laboratory frame of reference, the conservation equations written in divergence form are 

Dv 

mv+pl \ = 0, (1) 

mc 2 
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where v, D, m, E and p define, respectively, the fluid three-velocity, density, momentum density, total 
energy density and pressure. In the local rest frame the fluid can be described in terms of its Lorentz- 
invariant thermodynamic quantities: the proper rest mass density p, specific enthalpy h, and pressure p. 
The transformation between the two frames of reference is given by 

D = 1P, (2) 

m = P h l 2 -^, (3) 
E = ph 1 2 -p, (4) 

— 1/2 

where 7 = (l — t> 2 /c 2 ) is the Lorentz factor and c is the speed of light. The specific enthalpy, h, is 
related to the proper internal energy density, e, by 

h = — . (5) 

An equation of state (EoS) provides an additional relation between thermodynamic quantities and allows to 
close the system of conservation laws (1). Here we assume, without loss of generality, that the EoS expresses 
the specific enthalpy has a function of the pressure p and the specific proper volume r = 1/p, i.e. h = h(p, r). 
This allows us to define the sound speed c s as 

d P \ t 2 dh fdh V 1 

' -t , (6) 



de ) a h dr \dp 
where the derivative in equation (6) has to be taken at constant entropy s: 

dp 



dh\ s = [Tds + 

P 



= rdp. (7) 



We will consider only causal EoS, i.e., those for which c s < c. For such equations of state, the hyperbolic 
property of equations (1) is preserved (Anile 1989). In §3.4 we provide explicit expressions for the sound 
speed for different EoS suitable for relativistic flows. 

We describe the relativistic fluid in terms of the state vectors of conservative, U = (D, mi, m2, m^, E), 
and primitive, V = (p,v\,V2,Vs,p), variables. Here, and Vd (d = 1,2,3) are the projection of the three- 
momentum m and three- velocity v vectors along the coordinate axis, that is, = fid ■ m and Vd = fid • v. 

Equations (2)-(4) define the map U = U(V); the inverse relation gives V in terms of U: 

P=°, (8) 
V m (9) 



c 2 E + p ' 

p = Dh-y - E . (10) 

This inverse map is not trivial due to the non-linearity introduced by the Lorentz factor 7. Using equation 
(9) one can express the Lorentz factor as a function of the pressure: 

1 1-7^4- (ID 



7 2 (P) (E + pf 
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If D, m and E are given, equations (10) and (11) can be combined together to obtain an implicit expression 
for p: 

f(p) = Dh(p,T(p))i(p)-E-p = 0. (12) 

This equation must be solved numerically in order to recover the pressure from a set of conservative variables 
U. Notice that r = r(p) depends on the pressure p through t = -f(p)/D and that the specific enthalpy h is, 
in general, function of both p and r, h = h(p, r(p)). 

If the Newton-Raphson method is used as the root finder for pressure, the derivative of equation (12) 
with respect to p is needed. It can be expressed in terms of the derivatives of the specific enthalpy , already 
appearing in equation (6): 

#) =^-^f7? + MVl. (13) 



dp 'dp {E + pf \ ' 8t 

Equations (8)-(10), together with the pressure equation (12), have to be solved at each time step in order 
to provide the map V(U) needed in our numerical method. 



3. DESCRIPTION OF THE ALGORITHM 

Till now we considered the system of conservation laws (1]) written in divergence form. However, 
this differential representation is valid as long as the flow remains smooth. Since hyperbolic systems of 
conservation laws also admit discontinuous solutions, we should consider the more general integral form of 
the equations. We proceed as follows. 

We introduce a generic system of orthogonal curvilinear coordinates, (x 1 , x 2 , x 3 ), and define hd to be 
the unit vectors associated with the coordinate directions, d = 1, 2, 3. Orthonormality requires hd ■ fid' = 1 
when d = d' and hd ■ fid' = otherwise. The geometrical properties are then described in terms of the scale 
factors hi, hi and /13 which are, in general, functions of the coordinates, i.e., hd = hd{x 1 ,x 2 ,x 3 ). We divide 
our computational domain into control volumes, 

r x \ + L r x2 +i r x l+i 

AV ijk = 2 2 2 dV, dV ^hih 2 h 3 dx 1 dx 2 dx 3 , (14) 

J X 1 , J X 2 , J X 3 1 

*-3 i-h k -h 

where the index triplet k) corresponds to the hi, h 2 and h% directions. The integrals in equation (14) 
are evaluated between the lower ( , 1 ) and upper (x 1 , 1 ,x 2 , 1 , x, , 1 ) coordinate limits of the 

grid cell ijk, where 

x] +h -x)_ h =Ax], x) +h -x 2 _ h =Ax 2 , x 3 h+h -x\_ h =Ax\, (15) 

here Ax}, Ax 2 and Ax 3 , are the zone widths in the three coordinate directions. 

We restrict our attention to those coordinate systems for which the volume defined by equation (14) 
can be represented as the product of three different functions, each one depending on one coordinate only, 
i.e., 
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The three functions £ (x ), £ (x ), and £ (x ) are the generalized volume coordinates. These coordinates 
allow us to define the volumetric centroids £ 2 , ^ corresponding to the point (a;|,x 2 ,x|) as 

£ = e 1 ^ 1 ) = — ^ — ^ , (17) 



£ = e 2 (x|) = 2±S-^ , (18) 



2 

e k = e 3 (4) = • fc+2 2 • ^ . (19) 

Notice that in Cartesian coordinates, (^ 1 ,^ 2 ,^ 3 ) = (x,y,z), and the volumetric centroids coincide with 
the geometrical zone centers. In cylindrical geometry, (x , x 2 , x 3 ) — (r,(p,z), however, 

^M-y, C 2 (0) = 0, e(z) = z, (20) 
with the locations corresponding to the volumetric centroids given by 

/r 2 i+r 2 _i\ 5 7 -+i+0 7 -_i z k+ i+z k _i 

n = + 2 - 2 , & = J+2 - 3 2 , z fc - +2 * 2 . (21) 



Similarly, in spherical coordinates, (x 1 , x 2 , x 3 ) — (r,6,(j)), one obtains 

e(r) = j, e(0) = l-co S 9, C 3 (0)-0, (22) 

and the corresponding locations of the volumetric centroids are 

i 



- '< -7 \ _! /COS0 j+ i + COS0 ? _i \ _(j> k+ ^^ Vk _ 



r _ = f I , 9j = cos- 1 ( i^3_ £_L ] , = . ,23) 



The areas of the zone faces are expressed by A d — j n<i • dA, where dA is the area element dA — 
dA 1 ^ + dA 2 hi + dA 3 fi3, and dA d = dV / (ha dx d ). The result written in expanded form is 

Al ±hjk = H +i f k+h h 2 (xl ±i ,x 2 ,x 3 )h 3 (xl ±i ,x 2 ,x 3 )dx 2 dx 3 , (24) 

*/ 32 -i v 32 -i 



A\ j±1 _ k = / ' +2 f fc+2 h 1 {x 1 ,x 2 j± i,x 3 )h 3 {x 1 ,x 2 j± i,x 3 )dx 1 dx 3 . 

V 32 -y V 32 



(25) 



4 j;fe± i = / i ~ 2 1 2 2 h 1 (x 1 ,x 2 ,xl ± i)h 2 (x 1 ,x 2 ,xl ±h )dx 1 dx 2 . (26) 



Explicit expressions for the scale factors and face areas are given in Appendix A. 
Finally, the facial quadrature points located on the opposite bounding faces, 
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are used to compute the arc length separation, 

Al\ jk = J** +h h^x^xlxDdx 1 . (28) 

Similar expressions can be given for AlJ- k and Al\- k . 

The integral form of the equations can now be obtained by integrating the system (1) over a compu- 
tational volumes AVijk and over a discrete time interval At n = t n+1 — t n . Using Gauss's theorem, the 
evolutionary equations provide a relation between the volume averaged conserved quantities and the time- 
averaged surface integrals for the divergence terms: 

At 

'ijk J ijk 

At n 



D^k-m jk = -j^i> {Dv) n -dA, (29) 



fh n+1 — fh n — 

rn d,ijk " l d,i 3 k — 



ijk 



I (m d v) n -dA+ ( n d - (Vp + S G (mv)) n dV 

J ijk J ijk 



ijk -J ijk 

At 

ijk ~~ 111 ijk 



(30) 



K + k 1 -E? i k = -j={ r f (mr-dA, (31) 

LiVijk J ijk 



where d — 1,2,3 enumerate coordinate directions and the speed of light c has been set equal to unity, 
a convention to be used for the rest of this paper. Equations (29)-(31) express conservation of mass, 
momentum, and energy. In the equations above an overbar symbol denotes a volume average, 

5 «*=a^— / Dix^x^x^ndV, (32) 
while a bracket denotes a time-averaged quantity, 

i r tn+1 

(Dv) n =—J^ D(x\x 2 ,x 3 ,t)v(x\x 2 ,x 3 ,t)dt. (33) 

Notice that expressions (29)-(30) are exact and no approximation has been introduced so far. 

The vector Scimv) accounts for geometrical source terms arising from taking the divergence of the 
dyad mv in a curvilinear system of coordinates (explicit expressions for the geometrical source terms arc 
given in Appendix A). 



3.1. Evolution in Multi-dimensions 

Equations (29)-(31) are integrated by first solving the homogeneous problem, i.e. without the source 
term Sq in equation [30], separately from a source step. 

Let H At be the solution operator corresponding to the homogeneous part of the problem and S At be 
the operator describing the contribution of the source step over the time step At. Starting from the initial 
data, U n , the solution for the next two time levels is computed using operator splitting (Strang 1968): 

U n+1 = S At H At (U n ) , (34) 



whereas 



u n+2 = n At s At (cT +1 ) , 



(35) 
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Notice that this approach is second-order accurate in time, provided that the two operators have at 
least the same order of accuracy and the same time step is used for two consecutive levels. 



The homogeneous operator H. is based on the spatially unsplit fully corner-coupled method (CTU, 
Colella 1990; Saltzman 1994; Miller & Colella 2001, 2002). The conservative update results from acting H At 
on U n : 

n At (U n ) = U n + C 1 (t/£j , u£?) + C 2 ([/;+| , [/;+?) + £ 3 (u n k + + \ , U n k + _f) , (36) 

where J7™ = (L>yfc, m^fc, EV^)™ is a state vector of volume-averaged conserved quantities at time t = t n . 
Here we adopted the convention in which the integer- valued subscripts k are omitted when referring to 
three-dimensional quantities while the half integer values are used to denote zone edges. The same convention 
is used throughout the rest of the paper. 

The last three terms in equation (36), C d , represent one-dimensional flux-differencing operators, where 
d = 1,2,3 labels directions. For d = 1 (expressions for d = 2 and d = 3 are obtained by suitable permutations 
of indices) one has 



u: 



At 

AV 

At 

A! 1 



U. 



f 1 KJRi-; 



(37) 



where 



F\U) = 



is the flux vector, and 



P\U) 



DV! 

m\V\ 
m 2 v 1 
m 3 v 1 
mi 

/ \ 

P 





(38) 



(39) 



is the pressure term. The state vectors C^ i± i 2 are solutions to Ricmann problems with suitable time-centered 
left and right states properly computed at each facial quadrature point. For example, 



(40) 



where TZ(- 7 •) denotes the solution to the Riemann problem (see §3.3 for details). 



In the three-dimensional case, the conservative update (36) involves the following four steps. In the first 
predictor step we construct initial left and right states just as in the one-dimensional case. In the second 
predictor step the initial states are corrected for contributions from transverse directions to form secondary 
predictors. These secondary predictors are used in the third step to calculate the fully corner coupled states. 
Finally, the fluxes required by the conservative update follow from the solutions of Riemann problems with 
the fully corner-coupled states used as input. In what follows we present the details of this procedure. 
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The starting point for the construction of the input left and right states to the Riemann problems is the 
predictor step described in §3.2 below. In this step we use Taylor expansion to obtain edge- and time-centered 

- n+— ~ n+i 

estimates for the left and right states U i+ i L and U i+ i R . 



Next the following 6 secondary predictors are calculated: 



U 'h,s 




U k±±,S 


- U k±±,S 






U k±±,S 


- U k±±,S 


U 'ik,s 


- U i±i,S 




- tT h 



fr n+i 



Km 



u k+i 



.u 



n+i 



(41) 



(42) 



(43) 



where we adopt the convention that S — L at i + \ and S = R at i — i (similar formulae are obtained for 

The flux differencing operators on the right hand side of equations (41)-(43) are obtained by solving 
Riemann problems with the states obtained in the first corrector step, e.g., 



(44) 



In the third step we obtain the fully corner-coupled states, 

r n+i 1 



u: 

u 
u 



±i,S 
n+i 

fe±i,S 



U i±hs + 2 



c 3 (u' k 2 H ,u k 



'2 



u 



n+i 



1 



£ 2 {u'^ul 

;^ 3 Ki^^i)+^ o:+i^; 3 -i) 



u 



™+2 

fe±i,S ' 



(45) 

(46) 
(47) 



As in the previous step, quantities appearing in the flux differencing follow the solution of appropriate 
Riemann problems. For example, 



(48) 



The final step is the conservative update in which we use the fluxes obtained by solving Riemann 
problems with the fully corner-coupled states (eq. [45]-[47]). 

In three dimensions a total of 12 Riemann problems are required: 3 to obtain the secondary predictors 
(eq. [41]-[43]), 6 to obtain the fully corner-coupled states (eq. [45]-[47]), and 3 necessary for the conservative 
update (eq. [36]). The six secondary predictors, equations (41)-(43), are not required in two dimensions and 
the fully corner-coupled states become 



U. 



n + i 



T n +h 



+ 2 C 



1 „ 2 ( - n+i 

U i+i 



n+i 



(49) 
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^-^ + 5 £l K^)- <50 » 

Therefore, in 2D, the conservative update involves the solution of 4 Ricmann problems: 2 in the predictor 
step and 2 in the corrector step. 

Contribution from geometrical source terms is included by solving the ordinary differential equation 

d J = Sc(U), (51) 

with the solution from the previous step (either advection or possibly source term calculation if Strang 
splitting is used) used as the initial data. The solution is obtained with a second-order Runge-Kutta algo- 
rithm, and can be written in operator form as 

S A * (U) = \[U + U* + AtS G ([/*)] , 

2 (52) 
U* = U + AtS G (U) . 

Finally the choice of the time step At is based on the Courant-Friederichs-Lewy (CFL) condition 
(Courant et al. 1928): 

At = CaX ^{\^\^\^\)- ^ 
Here A^ ax is the fastest wave speed in the hd direction, and < C a < 1 is the limiting factor. 



3.2. The predictor step 

In what follows we provide a detailed description of the predictor step which constitutes the first step 
in the multi-stage procedure described above. For simplicity we describe the operator applied to the first 
coordinate direction; modifications required for other directions are straightforward. The predictor step is 
comprised of the reconstruction algorithm followed by the calculation of left and right states at the zone 
interfaces. These states are used as input data to the Riemann problems in the second predictor step (eq. 
[41]-[43]) and are required in the construction of the fully corner coupled states, equations (45)-(47) in 3-D 
and equations (49), (50) in 2-D. Throughout this section we make use of the following abbreviations: x = x 1 , 
Z = Z 1 ,h = h 1 and A = A 1 . 

During the reconstruction step a quartic polynomial is passed through the volume- averaged values £7". 
The polynomial provides interface values for the hydrodynamic variables, V" +i s (S = L,R) which arc 
fourth-order accurate in space in smooth regions of the flow. The interpolation may be done either in the 
volume coordinate £ (CW84) or in the spatial coordinate x (Blondin & Lufkin 1993). Here wc follow the 
former prescription and the details are given in Appendix B. 

Once the interpolated zone-edge values are provided, we define a one-dimensional piccc-wisc parabolic 
profile inside each cell i: 

for (54) 



v n (0 = vu iR + 



AV n + V£- I 1 - — 

61 I A& 
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whcrc 



V" , 



l - R ' 



r 61 



6 V? 



+ V U,n 



(55) 



Notice, however, that when the interpolation is done in the spatial coordinate x rather than in the 
volume coordinate £, equation (54) has to be modified by replacing £(x) with x. Moreover, the definition of 
Vei becomes geometry dependent (Blondin & Lufkin 1993). 

Second-order accuracy in time is obtained by considering the one-dimensional relativistic equations in 
quasi-linear form (derived for the one-dimensional case by Marti & Miiller 1996): 

dV 



Here 



V = 



f 


p 






Vl 






V2 






V3 






P 


) 



S(V) = 



vi d fldi 



A 2 <9£ \hdx 



( P \ 

v 2 c 2 s h 2 

V -phc 2 s j 

are, respectively, the vector of primitive variables and geometrical sources, while 

P 

Ml -^) 



(56) 







A <Y) = ^ 



o 



V 



V2C^ 
- 7 2 



phc 2 






Vl A 2 






hj 2 

2 



V 





Vl A 2 




phj 2 
ViV 2 (1 - 



phj 2 

vivs (1 ~ c 2 ) 
phj 2 

v 1 {l-c 2 s ) 



(57) 



is the Jacobian, A 2 = 1 - v 2 c 2 s and r/ 2 = 1 - vf - c 2 (v 2 + v%). 



n+4 



The left and right states V i+ i L and V i _ i R , are computed using the characteristic information carried 
to the zone edges £ i± i by each family of waves emanating from the zone center during a single time step 
At. Firstly, for each characteristic # we find the domain of dependence for each interface: 



where Ax # i and Ax^ i are defined through 



(58) 



x . 1 +Ai* 



/ 

</ X . , 1 — 



x. , i - Ax* , 



x 2 ,x\)dx = At max(0, —A*) , 
/i(x, x 2 ,x\)dx = At max(0, A*) . 



(59) 
(60) 
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Here # = {+, O^ 1 ' 2,3 - 1 , — } labels the characteristic field, and A # is the corresponding eigenvalue (Falle & 
Komissarov 1996; Donat et al. 1998), 



A~ = 



(l-c 2 s )-c s r) d,2,3) , 7«i (1 - c 2 ) + c s ?7 
' A =Bl ' = ^ ' 



(61) 



As in classical hydrodynamics, the characteristics consist of two genuinely non- linear fields =ff = +, — 
and three linearly degenerate fields # = O^ 1 ' 2 ' 3 -'. Also, note that in cylindrical or spherical coordinates, 
h = h}(x 2 ,x\) does not depend on x and therefore Ax^ ±1 = At max(0, ±\f )/h(x 2 , x\). 



Next, we calculate the average of V(£) over the domain of dependence lying to the left (for V i+ i L ) or 
to the right (for V i _i R ) of the interface: 



V* , 



i /■**+* 



(62) 



The zone averages are obtained by straightforward integration of the parabolic zone profiles (eq. [54]): 



i+i,L - V l +hL 



2A& 



£ i - f # 
AV n - Vr- I 1 - - * +hL 



V L, R = V U,n + 



2A& 



Avr + n i J i - - 



2?i-i,«-Si- 



A& 



(63) 



(64) 



Once the integrals are obtained, we use upwind limiting to select only those characteristics which 
contribute to the effective left and right states (Colella & Woodward 1984; Miller & Colella 2002). This 
requires expanding the matrix A(V) (eq. [57]) in terms of its eigenvalues and left and right eigenvectors. 
Specifically, upwind limiting for the left states consists in selecting those characteristics with positive speeds 
(i.e. emanating from the cell center towards the zone interface). Similarly, for the right states we consider 
only characteristics with negative speeds. The final result is: 



A*>0 



rf + = S< . 



At 



(65) 
(66) 



A*<0 



Here and r# are the left and right bi-orthonormal eigenvectors of A(V): 



( 



± 



\ 

PI 



2c s r] 



1 



2hc 2 . 



0(!. 2 .3) 



/ 1 \ 





1 



V 



he 2 



( 





1 



V2 



\ 



j 2 ph (1 — vf) 





V1V3 

o 
i 

V3 



j 2 ph (1 — vf) 



(67) 
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V 



1 

PI 

c s v 2 (jrjvi ± c s ) 
j 2 p(l~vf) 

c s v3 (rmvi ± c s ) 

j 2 p(l-v 2 ) 
he 2 



J 



1(1,2,3) 



/ 1 \ 





V J 





1 



v / 



/ o \ 




1 

v / 



(68) 



It can be verified that left and right eigenvectors are properly normalized, so that r # -Z # = 1 for # = #' 
and r# • =0 otherwise. Notice that, in the limit of vanishing transverse velocities, expressions (67) and 
(68) reduce to the one-dimensional left and right eigenvectors given by Marti & Miiller (1996) (apart from 
the normalization factor p7 2 /c s int and I ). 



3.3. Solution of the Riemann Problem 

The Riemann solver describes the evolution of a discontinuity separating two arbitrary constant hydro- 
dynamic states, V l and V r: 

V(x,t = 0) = \ Z L !° r X< 1> (69) 

\ V R for i>0. K ' 

Like in the Newtonian case, the solution is self-similar (i.e. it is a function of x/t alone), and the decay 
of the initial discontinuity gives rise to a three-wave pattern (see Fig. 1): a contact discontinuity separating 
two non-linear waves, either shocks or rarefactions. The contact discontinuity moves at the fluid speed and 
both the normal velocity v\ and pressure p are continuous across it, while density p and tangential velocities 
i>2,3 generally experience jumps. Across a shock wave, however, all the components of V can be discontinuous 
and their values ahead and behind the shock are related by the Rankine-Hugoniot jump conditions. Inside a 
rarefaction wave, on the contrary, density, velocity and pressure have smooth profiles given by a self-similar 
solution of the flow equations. In the multidimensional case, the relation between the hydrodynamical states 
at the head and the tail of the rarefaction wave can be cast in one non-linear ordinary differential equation 
(Pons et al. 2000) the solution of which is usually computationally expensive. Therefore our solver relies 
on the "two-shock" approximation: a Riemann problem is solved at each zone interface by approximating 
rarefaction waves as shocks. This approach has been extensively used in Godunov-type codes (e.g., Colclla 
& Woodward 1984; Woodward & Colella 1984; Dai & Woodward 1997). 

The jump conditions for the left (S — L) and right-going shock (S = R) can be written as (Taub 1948; 
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Pons et al. 2000): 

' ' = -CM , (70) 



[hiv 1 ] = {[p], (71) 
[h 1V2 ]=0, (72) 
[hjv 3 }=0, (73) 

hi - ^] = C[pvi] , (74) 

where [q] = q — qs denotes the difference between the pre-shock (qs, S = L,R) and the post-shock state q, 
C = lsh/j, Ish = (1 — u sfc) -1 ^ 2 is the Lorentz factor of the shock, v s h is the shock velocity and j is the mass 
flux across the shock. 

We note that in relativistic flows the tangential velocities can exhibit jumps across the shock (Pons 
et al. 2000; Rezzolla et al. 2003). This is one of the major differences from the classical case caused on 
one hand by the coupling between different velocity components through the Lorentz factor (kincmatically 
relativistic regime where the flow speeds become comparable to the speed of light) and on the other hand 
by the specific enthalpy being tied to the Lorentz factor (thermodynamically relativistic regime when the 
pressure p is comparable to or greater than the rest mass energy density p) . Both characteristics are absent 
from Newtonian hydrodynamics. 

The solution of the Riemann problem is obtained as follows. First, we combine equations (70)-(74) to 
obtain the Taub adiabat (Taub 1948), the relativistic version of the Hugoniot shock adiabat. The adiabat 
relates thermodynamic quantities - enthalpy h, pressure p, and specific proper volume r - ahead and behind 
the shock: 

[h 2 ] = (hr + h S T S ) [p] . (75) 

The above equation together with an EoS, h = h(p,r), is solved to obtain r and h as functions of the 
post-shock pressure p (explicit solutions for selected equations of state are provided in §3.4). These solutions 
together with the definition of the mass flux, 

j 2 = y 2 ah D 2 s (v ah -vs) 2 , (76) 
lead to a quadratic equation for £. That equation can be solved to obtain 

V S V S ± y/Vi + {l^^ s )/p{p,S) 

where V$ — 1/Ds, and for the mass flux we use the expression (Marti & Miiller 1994, 1996; Pons et al. 2000) 

> 2 =~wv <78) 

with h and r being functions of p and S. In writing equation (77) and in what follows we adopted the 
convention of Pons et al. (2000) with the plus (minus) sign corresponding to S = R (S = L). Explicit 
expressions for j 2 (p, S) for four different EoS are given in §3.4. 



C± (P, S) = v * 1 \ b " J ^ ' , (77) 

1 Vc 



-14- 



We now express the post-shock velocity as a function of the post-shock pressure. We use equations (77) 
and (71) with hj in the post-shock state given by equation (74) and 1/D given by equation (70) (Marti & 
Miiller 1996): 

v(p, s) = h s-ysvs + \pK(P,s) (79) 



hsls + \p] (vsC(p, S) + V s S j 



The solution to the Riemann problem is now completely parameterized in terms of the post-shock 
pressure p. The solution for the normal velocity and pressure, v* and p*, has to be continuous across the 
contact discontinuity. We achieve this by imposing the condition v(p, R) = v(p, L) as shown in Figure 2. 
Since this equation cannot be solved in closed form, an iterative scheme has to be employed. Here we propose 
a Newton- Raphson algorithm. Given the (n)-th iteration to p* , the (n+1) approximation is found through 

(n+l) = (n) _ <P (n \L) - v(pW,R) 

P P v'(pM,L)-v'(p^,R)- [ 1 

The explicit expression for the derivative v'(p( n \S) — dv(p,S)/dp\ p — p ( n ) is found by differentiating 
equation (79) with respect to the post-shock pressure p: 

dv(p S) (C(P. S) + W(P, S)) (l - v s v(p, S)) - v{p, S)V S 
v'( P ,S) = 7 ; = ^ L . (81) 



dp h S 7s + \p}(CiP,S)vs + Vs) 



Here we take advantage of the fact that the derivative of £'(p,S) = d£(p,S)/dp, appears only through 
[p]C(p,S). This term can be eliminated by differentiating equation (77) with respect to the post-shock 
pressure and multiplying the result by [p] . With some algebra we obtain 

1 d(hT)/dp+l/j 2 ( P ,S) 

2 C(p,S)(l-v 2 s )-V s v s 

Analytic expressions for d{hr)/dp for selected EoS are given in §3.4. When no analytic expression is available 
or is expensive to compute, we replace the explicit expression for v'(p, S) with a numerical approximation 
based on the most recent iteration: 



rWKS)*^- 

To start the pressure iteration we use 



(p^, S) - vip^-V, S) 



p(n) — p(n-l) 



„, fes) = ± Vl±t3SW), (84) 
hsls 

with the + (— ) sign for S = R (S = L). Note that the term j 2 (ps,S) appearing in the above equation 
cannot be computed using equation (78) since this expression becomes ill-defined when p — > ps- We address 
this problem in §3.4. 

The iteration process stops once the relative change in pressure falls below 10~ 4 -10~ 6 . Usually no more 
than 6-7 iterations are needed to achieve convergence even for the most extreme test cases presented in §4. 

Once v* and p* have been obtained, the tangential velocities v^g and v^ s on both sides of the contact 
discontinuity are computed using equations (72) and (73) with hj in the post-shock state given by equation 
(74) and 1/D given by equation (70): 

* JshsV2S,3S ,oc\ 

U2S ' 3S " lsh s + (P* - Ps) (v s C(P*,S) + V S ) • 1 ' 
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Finally we use equation (70) to obtain the proper densities p* s . 

The solution to the Riemann problem is now represented by the four regions in the (£, t) plane: Vl, V* l , 
V* R and V a (Fig. 1). The numerical fluxes required in §3.1 are functions of the state vector 1Z(Vl, Vr), 
computed by sampling the solution to the Riemann problem at the £/t = axis. To this purpose, we define 
c = — sgn(v*); it follows that 



n(v L ,v R ) = { 



v s 
v s 

A s -A s 



where we set 



V* s ,Vs 



VI, V L 
V* R ,V R 



if aX* s > , 
if a\ s < , 

if a\* s < < a\ s ■ 

if a < , 
if a > . 



(86) 



(87) 



The last condition in equation (86) holds when the solution is sampled inside a rarefaction fan. In this 
case we obtain TZ(Vl, V r) by linearly interpolating the fan between the head and the tail. 

We compute X* s and As according to whether the wave separating Vs from V s is a shock or a rarefaction. 
If p* > ps then the wave is a shock, and we set both As and X* s to be equal to the shock speed 



As = Ac 



1 



+ v s ■ 



(88) 



S D s ((p*,S) 

If, on the other hand, p* < ps the wave separating Vs from V* s is a rarefaction. In this case we first define 



Ac 



X~(V* S ) if a < 0, 

{ X+(V* R ) if a>0, 

and then enforce the condition A^ < A£ (if a < 0) or X R > X* R (if a > 0) by setting 

(x-(V L ),Xl) if a<0, 



(89) 



As 



^ max ^A+ (V R ) , X* R 



if cr > . 



(90) 



The expressions for the maximum and minimum wavespeeds A + and A are given by equation (61). 



3.4. Equation of State 

The results from the previous section have shown that our Riemann solver can be written in terms of 
general expressions involving the mass flux (eq. [78]), the solution to the Taub adiabat (eq. [75]), and the 
derivative d(hT)/dp (eq. [82]). Their explicit forms, however, depend on the particular choice of the EoS. In 
what follows we consider four different equations of state that can be cast in the form f(h, 6) = 0, where 
6 = pr is a temperature-like variable. EoS properties are conveniently described in terms of the function 

r.. ( e, s MS_l. (91) 
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Finally, we also derive some useful expressions in constructing our special relativistic hydrodynamics code. 

The first equation of state is represented by the ideal gas EoS (hereafter labelled ID) characterized by 
a constant polytropic index T. In this case 

/i(G) = l + f ^ T 6. (92) 

For this EoS, T r is simply a constant and varies between 5/2 (for T — 5/3) and 4 (for T — 4/3). The 
ideal gas EoS is very popular in non-relativistic fluid dynamics and has been largely used also in studying 
relativistic flows. However, equation (92) is not consistent with relativistic formulation of the kinetic theory 
of gases and admits supcrluminal wave propagation when T > 2 (Taub 1948). In his fundamental inequality, 
Taub proved that the choice of the EoS for a relativistic gas is not arbitrary. The admissible region is defined 

by 

(h-e)(h-4Q) > 1, (93) 

shown in the upper-left portion in Figure 3. Note that for a constant T-law gas, the above condition is always 
fulfilled for T < 4/3 while it cannot be satisfied for T > 5/3 for any O > 0. Furthermore, by inspecting 
the expressions given in Table 1 one can see that the ideal gas EoS admits supcrluminal sound speed when 
T > 2 (corresponding to T r < 2). The exact form of relativistic ideal gas EoS was given by Synge (1957). 
For a single specie the relativistic perfect {RP) gas is described by 

h K3il/e) (94) 

where K 2 and K3 are the modified Bessel functions of the second and third kind, respectively. For the RP 
EoS T r (9) reduces to 5/2 in the limit of low temperature (6 — > 0) while for an extremely hot gas (6 — ► 00) 

r r (e)^4. 

Unfortunately no simple analytical expression is available for this EoS and the thermodynamics of 
the fluid is entirely expressed in terms of the modified Bessel functions (Falle & Komissarov 1996). The 
consistency comes at the price of extra computational cost since the EoS is frequently used in the process of 
obtaining numerical solution. Below we offer two more formulations which use analytic expressions and are 
therefore more suitable for numerical computations. 

Sokolov et al. (2001) proposed a simplified EoS which has the correct asymptotic behavior in the limit 
of ultrarelativistic temperatures (6 — ► 00). Their "interpolated" EoS {IP) 

h(h - 46) = 1 , (95) 

has the attractive feature of being simple to implement and is more realistic than the ideal EoS (eq. [92]). 
However, the IP EoS differs from the exact RP EoS in the limit of low temperatures where T r (9) — > 2 
rather than T r (6) — > 5/2. Besides, the IP EoS does not satisfy Taub's inequality (eq. [93]) as can be seen 
from Figure (3). 

Here we propose a new EoS which is consistent with Taub's inequality for all temperatures and it has 
the correct limiting values. Our newly proposed EoS (TM) is obtained by taking the equal sign in equation 
(93) to obtain 

(h-G)(h-4G) = 1. (96) 



The TM EoS is quadratic in h and can be solved analytically at almost no computational cost. We 
point out that only the solution with the plus sign is physically admissible, since it satisfies ft, > 1 for all 9 
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and it has the right asymptotic values for T r (9). Direct evaluation of r,.(0) shows that the TM EoS differs 
by less than 4% from the theoretical value given by the relativistic perfect gas EoS (eq. [94]). 

For a given EoS, our Ricmann solver requires the solution of the Taub adiabat (eq. [75]), followed by 
the calculation of the mass flux (eq. [78]) and the derivative of hr with respect to pressure appearing in 
equation (82). 

For the ID EoS (eq. [92]) and the TM EoS (eq. [96]) the Taub adiabat reduces to a quadratic equation 
in x = h and x — [hr], respectively: 

ax 2 + bx + c = . (97) 
For the ID EoS the coefficients a, b and c are given by 

= 1-=^-, b=^, c = -h s (hs + Ts\p]), (98) 
i r p i r p 

while for the TM EoS one has 

a = p S (3p + Ps), b= -h 2 s (3p + 2p s ) - h S T S (3p 2 - 7pp s - 4p|) - [p] , 
c = -h S T S \p] (h 2 s + 2h S T SP + 2) . 



(99) 



The solution is even simpler for the IP EoS (eq. [95]) but has to be found numerically for the RP EoS (eq. 
[94]), see Table 1. 

Straightforward evaluation of the mass flux, equation (78), leads to numerical difficulties whenever 
[p] — > because at the same time also [hr] — > 0. For all the equations of state presented here, with the 
exception of the RP EoS, one can find alternative expressions that do not suffer from this drawback. The 
procedure is to factor out [p] from [hr]. The results, obtained after some algebra, are given in Table 2. When 
such a factorization is not possible, the limiting value of j 2 can be found by noticing that 



(100) 



^5)=-W/*l b ]=o- 
Since from the Taub adiabat one also has dh/dp\[ p ] =0 = ts , we obtain the final result 

j 2 ( Ps ,S) ee lim f(p,S) = * ksPs , (101) 

[p]-"0 T s h s psTs + hs(l - h s ) 

where hs — dh(Q)/dQ\^ =Q . Finally, the derivative d(hr)/dp is obtained by straightforward differentiation. 



4. NUMERICAL TESTS 

We use several problems to verify the correctness of the algorithm and its implementation. We consider 
one- and multi-dimensional setups and problems close to real astrophysical applications. 

The first two problems are of particular interest in the context of Ricmann problem as they involve the 
decay of an initial discontinuity into three elementary waves: a shock, a contact and a rarefaction wave. 
In this case an analytical solution can be found (see Pons et al. 2000; Marti & Miiller 1994, and references 
therein) . Shock-tube problems are basic verification tests useful in assessing the ability of the Riemann solver 
to correctly describe evolution of simple waves. Thanks to their relative simplicity, one dimensional problems 
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are also excellent benchmarks that can be used for comparison of different algorithms and implementations 
(Calder et al. 2002). 

We also consider one of the two-dimensional Riemann problems originally introduced by Schulz et al. 
(1993) and later on adopted for relativistic flows by Del Zanna & Bucciantini (2002). The two-dimensional 
Riemann test problem is relatively easy to set up and is characterized by increased complexity. Four ele- 
mentary waves, two shocks and two contact discontinuity, are present in the initial conditions; their mutual 
interaction is also calculated. 

We propose a three-dimensional test involving the reflection of a spherical relativistic shock (Aloy et al. 
1999) in cartesian coordinates. 

Finally, we consider two astrophysical applications: pressure-matched light relativistic jet and evolution 
of the jet in stratified medium. The latter problem can be relevant in the context of gamma-ray bursts and 
their afterglows. 

4.1. Relativistic Shock Tubes 

In one-dimensional shock-tubes a diaphragm initially separates two fluids characterized by different 
hydrodynamic states. At t = the diaphragm is removed and gases are allowed to interact. We consider 
two cases and assume an ideal equation of state with T = 5/3. The calculations are performed in cartesian 
geometry. The domain covers the region < x < 1 and the discontinuity separating the two states is initially 
placed at x — 0.5. Free outflow is allowed at the grid boundaries. 

In the first problem (Pl) density and pressure to the left of the interface are pl — 10 and pl = 40/3. 
To the right of the interface we have pn — 1 and p^ = 2/3 • 10~ 6 . The gas is initially at rest. This problem 
has been extensively used by many authors (see Marti & Miiller 2003, and reference therein); our results at 
t = 0.36 are shown in Figure 4. 

In this problem a rarefaction wave develops and moves to the left while a contact discontinuity separating 
the two gases and a shock wave are both propagating to the right. The shocked gas accumulates in a thin 
shell and is compressed by a factor of w 5. 1 Relativistic effects are mainly thermodynamical in nature. The 
enthalpy of the left state is appreciably greater than in the non-relativistic limit, h — 10/3 > 1, while the 
Lorentz factor is rather small (7 ~ 1.4). With 400 equidistant zones and C a = 0.9 our algorithm is able to 
reproduce the correct wave structure with considerably high accuracy, as shown in Figure 4. The shock and 
the contact discontinuity are captured correctly. The contact is smeared only over 2-3 two zones while shock 
profile occupies 3-4 zones. The thin dense shell is well resolved. The rarefaction fan is also well reproduced 
but a slight overshoot in velocity and undershoot in density appears at the tail of the rarefaction. A direct 
comparison with the analytical solution (see Table 3) shows that, at the resolution of 400 zones, the error 
measured in the L — 1 norm is smaller than 5%. We obtained first order convergence with increasing grid 
resolution, as expected for discontinuous problems. 

In the second one-dimensional test (P2) the gas to the left of the discontinuity has pressure pl = 10 3 
while pr = 10~ 2 . The proper gas density and velocity are uniform everywhere and equal to p = 1 and 
v x = 0, respectively. The initial conditions are therefore similar to those of Pl, but the internal energy 
of the left state is greater than the rest-mass energy, h ~ 2.5 • 10 3 » 1, resulting in a thermodynamically 



1 In relativistic hydrodynamics, the density jump across a shock wave can attain arbitrary large value. 
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relativistic configuration. We also include the effects of non vanishing tangential velocity. The test is run 
for nine possible combinations of pairs v y x and v v ^r constructed from the set of {0, 0.9, 0.99} (Pons et al. 
2000). The results at t = 0.4 are shown in Figure 5. 

The decay of the initial discontinuity in P2 leads to the same wave pattern as in PI, but the changing 
tangential velocity has a strong impact on the solution. In absence of shear velocities (Fig. 5a), the high 
density shell is extremely thin posing a serious challenge for any numerical scheme. With 400 zones and 
C a = 0.4 the post-shock density is underestimated by ~ 25%. The error is reduced to ~ 6% when the 
resolution is doubled. 

The solution changes drastically when tangential velocities are included due to the increased steepness 
of the shock adiabats, dp/dv. This can be seen from equations (79) and (81) setting vs = 0. This is a purely 
relativistic effect resulting from the coupling of the specific enthalpy (which is independent of transverse 
velocities) and the Lorentz factor in the upstream region. Specifically, a higher tangential velocity in the 
right state has the effect of increasing the post-shock pressure, so that, even if the effective inertia of the 
right state is larger (due to the higher Lorentz factor), the two effects compensate and the shock speed is 
only weakly dependent on v y R. As a net result the mass flux increases with the Lorentz factor of the right 
state 7 fl . This leads to the formation of a denser and thicker shell (cases b, c, and to smaller degree e and / 
in Fig. 5). 

The opposite behavior is displayed when the tangential velocity is large in the left state (cases d and g), 
leading to a lower post-shock pressure and smaller normal velocity. As a consequence, the shock strength 
and the mass flux through the shock are largely reduced. 

Figure 5 presents the nine numerical tests together with the corresponding exact analytical solutions 
(Pons et al. 2000). One can clearly notice an excessive smearing of the contact discontinuity with increasing 
v v l. The situation does not improve even when the exact Riemann solver is used and the evolution of 
rarefaction waves is calculated with greater accuracy. We speculate that the reason for excessive smearing of 
the contact discontinuity is that our numerical scheme is not able to properly resolve waves right after breakup 
of the initial discontinuity. We expect that other currently existing shock-capturing schemes experience 
similar diffuculties in this case, but such a detailed study is beyond the scope of the present publication. 

4.2. Two-dimensional Riemann Problem 

Two-dimensional Riemann problems have been originally proposed by Schulz et al. (1993) in the context 
of classical hydrodynamics (Lax & Liu 1998). This set of problems involves interactions of four elementary 
waves (shocks, rarefactions, and contact discontinuities) initially separating four constant states. A relativis- 
tic extension of a configuration involving two shocks and two contact discontinuities was recently presented 
by Del Zanna & Bucciantini (2002). Here we consider this test but with slightly modified initial conditions 
that reproduce a simple wave structure (the conditions adopted by Del Zanna & Bucciantini (2002) do not 
yield elementary waves at every interface). 

We define four constant states in the four quadrants: x,y > (region 1), x < 0,y > (region 2), 
x,y < (region 3), and x > 0,y < (region 4). The hydrodynamic state for each region is given in the 
left panel in Figure 6 with p\ = 5.477875 • 10~ 3 and p\ ~ 2.762987 • 10~ 3 The computational domain is the 
box [—1, 1] x [—1, 1]. We used an ideal EoS with T = 5/3. The integration is carried out with C a = 0.4 till 
t = 0.8. 
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The initial conditions result in two equal-strength shock waves propagating into region 1 from regions 
2 and 4 (right panel in Fig. 6). The two contact discontinuities bound region 3. By t — 0.8 the unperturbed 
solution is still present in part of region 1 and 3. Most of the grid interior is filled with the shocked gas. 
The two shocks continuously collide in region 1. Curved transmitted shock fronts bound a double-shocked 
drop-like shaped region located along the main diagonal. At the final time the symmetry of the problem is 
well preserved with a relative accuracy of about 10 -9 in density. We also note that our nonlinear Riemann 
solver can capture stationary contact discontinuities exactly, without producing the typical smearing of more 
diffusive algorithms (Del Zanna & Bucciantini 2002). 



4.3. Relativistic Spherical Shock Reflection 

The initial configuration for this test problem consists in a cold (p = 0), uniform (p = 1) medium with 
constant spherical inflow velocity v in and Lorentz factor -f in . At t = the gas collides at the center of the 
computational domain and a strong shock wave is formed. For t > the shock propagates upstream and the 
solution has an analytic form given by (Marti et al. 1997): 

for r > v s t 

(102) 

for r < v s t 
where 

r + 1 r 

is the compression ratio and 

r - 1 

V s = —rlin\v in \ (104) 

7 + 1 

is the shock velocity. Here a — 0, 1, 2 for cartesian, cylindrical and spherical geometry, respectively. Behind 
the shock wave (r < v s t), the gas is at rest (i.e. v = 0) and the pressure has the constant value p(r, t)(j in — 
l)(r — 1). Conversely, in front of the shock all of the energy is kinetic and thus p — 0, v = v in . 

This test problem has been proposed by Aloy et al. (1999) in three-dimensional cartesian coordinates 
to test the ability of the algorithm in keeping the spherically symmetric character of the solution. For 
numerical reasons, pressure has been initialized to a small finite value, p = 2.29 • 10~ 5 (r — 1), with T = 4/3. 
The computational domain is the cube [— 1, l] 3 and the final integration time is t = 2. 

Figure 7 shows the results for i> in = —0.9 on 101 3 computational zones; the Courant number is C a = 0.4. 
The relative global errors (computed as in Aloy et al. 1999) for density, velocity and pressure are, respectively, 
10.7%, 0.7%, 16.4%. Numerical symmetries (computed as max(p(xi, yj, 0) — p(0, yk, Zj j) are retained to 10~ 7 
for this case. From the same Figure we notice that our solver suffers from the wall-heating phenomenon, 
typical of such schemes. 

We also investigated the ultra-relativistic regime by increasing the inflow velocity according to v in = 
v—1, where v = 10 _1 , 10~ 3 , 10~ 5 , 10~ 7 , the latter corresponding to a Lorentz factor of j in w 2236. Following 
Aloy et al. (1999) we perform this set of simulations on 81 3 zones. The relative global errors are given in 
Table 4. 



P(r,t) 



1 + \Vi- 



1 + 
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4.4. Multi-dimensional Applications 

Relativistic hydrodynamics is expected to play an important role in many phenomena typical in high 
energy astrophysics. In this Section we present the results of application of our new algorithm to the problem 
of propagation of relativistic jets. In the first problem we consider a highly supersonic pressure-matched 
light jet propagating in a constant density medium. In our second example we study the evolution of a light, 
high-Lorentz factor jet in a stratified medium, a situation frequently considered in the context of cosmological 
gamma-ray bursts. 

4- 4-1- Axisymmetric Jet 

We consider the propagation of an axisymmetric relativistic jet in 2-D cylindrical coordinates. The 
computational domain covers the region < r < 15 and < z < 45 jet radii. At t — the jet is placed on 
the grid inside the region z < 1, r < 1 with velocity v z = 0.995 parallel to the symmetry axis and jet beam 
density pi, = 1. The ambient medium has density p m — 100 and pressure is uniform everywhere, p — 0.005. 

We use the TM equation of state, equation (96), which for the current conditions yields a relativistic 
beam Mach number of jVz/ (c s -f Cs ) 110.5 (corresponding to a classical Mach number Vb/c s s=a 11). The 
computational domain is covered by 360 x 960 equally spaced zones in the n r and h z direction respectively. 
At the symmetry axis we apply reflecting boundary conditions. Free outflow is allowed for at all other 
boundaries with the exception of the jet inlet (z = 0, r < 1) where (p, v,p) are kept at their initial values. 
The integration is performed with C' a = 0.4. 

We follow the evolution until t = 80 (in units of jet radius crossing time at the speed of light). By 
that time the jet has developed a rich and complex structure (Fig. 9). The morphology of the jet cocoon 
is determined by the interaction of the supersonic collimated beam (inside r < 1) with the shocked ambient 
gas. The beam is decelerated by the Mach disk and is terminated by the contact discontinuity where the 
beam material is deflected sideways feeding the cocoon. The interface between the jet backflowing material 
and the ambient shocked gas is Kelvin-Helmholtz unstable resulting in the formation of vortices and mixing. 
For the current conditions the jet cocoon shows a mix of turbulent regions close to the beam and relatively 
smooth outer sections traversed by a number of weak sound waves. The beam shows 3-4 internal shocks in 
its middle section and another strong internal shock about 10 jet radii behind the jet head. 

Our results favorably compare to the previous study of Marti et al. (1997). Although in our simulation 
we do not use an ideal EoS, the equivalent adiabatic index, defined through equation (92) as T eq = (h — 
l)/(/i — 1 — O), always remains very close to 5/3 making direct comparison with other studies possible. 
Specifically, our choice of parameters places the presented model between models C2 and C3 of Marti et al. 
(1997). Figure 8 shows the position of the jet head as a function of time. We find that during the early 
evolution (t < 40) the jet head velocity agrees very well with the one-dimensional theoretical estimate given 
by Marti et al. (1997). At later times, 40 < t < 60, the jet enters a deceleration stage after which (t > 60) 
the velocity becomes again comparable with the theoretical prediction. 

The propagation efficiency, defined as the ratio between the average head velocity and its one-dimensional 
estimate, is found to be w 0.92. This result is consistent with the previous results obtained for highly- 
supersonic light jets with large adiabatic index for which beam velocities tend to have propagation efficiency 
close to 1 (see Table 2 in Marti et al. 1997). 
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4-4-2- Gamma Ray Burst 

We propose a simplified model for the propagation of a relativistic jet through a collapsing, non-rotating 
massive star. Such objects are believed to be sources of long-duration gamma ray bursts observed at cosmo- 
logical distances (Frail et al. 2001; Woosley & MacFayden 1999; Aloy et al. 2000). Our choice of parameters 
closely reflects that of model JB of Zhang et al. (2003). 

We adopt spherical coordinates (r, 9). The domain has 320 uniformly spaced zones between 1 < r/r < 
11 and 500 geometrically stretched zones in the range 11 < r/r < 500. Here r = 2 • 10 s cm is the inner 
grid boundary. The grid in the 9 direction consists of 180 uniformly-spaced zones for 0° < 9 < 30°, while it 
is geometrically stretched with 100 zones covering the region 30° < 9 < 90°. 

For the sake of simplicity, density and pressure distributions in the stellar progenitor are given by single 
power-laws, 

Prn = PO (J^j , Pm=Pa (J^j , (105) 

while the velocity is zero everywhere. We adopted a = 1.77, b = 2.6, po = 10 6 gr/cm 3 , and po = 8.33 • 10 23 
dyn/cm 3 to match the progenitor model of Zhang et al. (2003). To account for gravity we included an 
external force corresponding to a 5M Q gravitational point mass placed at r = 0. This is accounted for 
during the source step of our algorithm. 

The jet is injected at r — r , < 9 < 9 where 9 is the jet half-opening angle. The jet is characterized 
by its Lorentz factor jb, energy deposition rate (jet power) 

E b = Pbjb(hjb - l)v b 2nrl(l - cos6» ) , (106) 
and kinetic to total energy density ratio 

= Pblb (n b i) (107) 

Pblb{hb^b - 1) ~Pb ' 

where pb, hb and pb are the proper density, specific enthalpy and pressure of the beam, respectively. Notice 
that the total energy density considered here does not include the rest mass energy. Equations (106) and 
(107) are used to express pressure and density as functions of the jet power, fractional kinetic energy density, 
Lorentz factor, and jet half-opening angle. For the present application we use Eb = 10 51 erg s _1 , / = 0.33, 
7h = 50 and 9o = 5°. The TM equation of state is adopted and the CFL number is C a = 0.6. Reflecting 
boundary conditions are used at 9 — and 9 = 90°, while we allow for a free outflow at the outer grid 
boundary (r/ro = 500) and at the inner boundary outside the injection region. We follow the evolution until 
the jet has reached the outer boundary at r = 10 11 cm. 

The results of our simulation are shown in Figure 1 1 . During the early stages (t < 1 s, left panel in Fig. 
11) the jet beam is quickly collimated by the high pressure of the near stellar environment, p m /Pb ~ 263.5. 
By t = 3.5 s (right panel in Fig. 11), the jet has cleared a low-density, high-velocity thin funnel along the 
polar axis and remains narrowly collimated with a very thin featureless cocoon. The latter property reflects 
the fact that ultra-relativistic jets arc less prone to Kelvin- Helmholtz instabilities (Ferrari et al. 1978; Marti 
et al. 1997). This behavior is typical for supersonic light jets characterized by low Mach numbers (Mb ~ 1.82 
and Pb/ p m \r=r ~ 1-85 • 10~ 5 for the present case). One should also be aware that our simulation includes 
a certain amount of numerical diffusion which limits the resolution of small scale structures in our model. 
This is specially true at large radii where the resolution of our mesh is coarser (the zone width increases 
geometrically with radius) and the geometry is diverging (spatial resolution decreases laterally) . 
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Our results stay qualitative in agreement with the results of Zhang et al. (2003) with only minor 
differences. Specifically, the jet in our model propagates about 50% faster for t < 0.7 s but reaches a similar 
distance by the final time (t = 3.5 s). The morphological differences between the two models are small, with 
the current model having a slightly narrower beam and a less prominent cocoon. These discrepancies likely 
result from a combination of factors including differences of the numerical methods, initial conditions, and 
physics (we used simplified equation of state and neglect radiation effects). Overall, however, the evolution 
of the jet in both models is very similar. 

5. CONCLUSIONS 

We presented a high-resolution numerical scheme for special relativistic multidimensional hydrodynam- 
ics in general curvilinear coordinates. A finite volume, Godunov-type formulation is used, where volume 
averaged conserved quantities are evolved in time by solving Ricmann problems at each time step. The 
solver takes into account non- vanishing tangential velocities at each zone interface and assumes that the two 
non-linear waves are shocks (i.e. "two-shock" approximation is used). This greatly reduces the computa- 
tional cost and turns out to provide a reasonable approximation in the limit of weak rarefaction waves (as 
long as the time integration is done explicitly so that the time step is relatively small). The solution to the 
Ricmann problem is found iteratively and a new method of incorporating a general EoS is presented. 

We considered four different equations of state suitable for relativistic hydrodynamics. A novel simple 
analytic formulation for the relativistic perfect gas EoS has been presented. Our new equation of state 
recovers the exact solution (Synge 1957) with accuracy better than 4%. This formulation is consistent with 
a special relativistic formulation of the kinetic theory of gases and shows the correct asymptotic behavior in 
the limit of very high and very low temperatures. Since our EoS is given by a simple analytic expression, 
the computational cost of the solution of the Riemann problem is significantly reduced. 

Multidimensional integration is done with the fully coupled corner-transport upwind method (Colella 
1990; Saltzman 1994; Miller & Colella 2002). Preserving symmetries of the problem is often difficult especially 
when directionally split advection algorithm is used and may require applying special procedures (Aloy et al. 
1999). Our choice of an unsplit integrator makes the presented method free of such problems and symmetries 
of the flow are perfectly preserved. 

The calculation of the numerical fluxes requires solving 4 Riemann problems in two dimensions and 
12 Riemann problems in three dimensions per cell per time step. This compares to 6 Riemann problems 
to be solved for the unsplit second-order Runge-Kutta schemes (in three dimensions) which in practice, 
however, usually require smaller time steps. Second-order accuracy in time is achieved by using characteristic 
projection operators. 

Our implementation is verified in the case of strongly relativistic flows with Lorentz factor in excess of 
2 • 10 3 . Purely one-dimensional test problems with non- vanishing tangential velocity are used to demonstrate 
the formal correctness of the method. The implementation is verified in two dimensions using a Riemann 
problem with two shocks and two contact discontinuities. Performance of the new method is illustrated 
using two astrophysical applications: a pressure-matched light relativistic jet with Mach number « 11 and 
r m 10 in 2-D cylindrical geometry and the propagation of a highly relativistic (r > 50) jet through the 
stellar atmosphere in spherical geometry. The latter problem is motivated by studies of relativistic jets in 
collapsars (Zhang et al. 2003). 
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A. ORTHOGONAL CURVILINEAR COORDINATES 

In what follows we show how the geometrical source term Sg (eq. [30] ) can be calculated for an arbitrary 
orthogonal system of coordinates. We also provide explicit expressions for the scale factors, face areas and 
zone volumes as required by the conservative finite volume formulation presented in §3. Expressions are 
given for the most commonly used coordinate systems: cartesian, cylindrical and spherical. 

Consider the rank-two tensor, T = J2 a b Tabn a nb- The divergence of T is the vector 

where Sg(T) is the source term contributed by those versors that do not have fixed orientation in space. 
The source term vector can be expressed as 

c (r\ ST ( Tld dfld _l T2d dfld _l T3d dft A (A o\ 



d=i 

with components 



a « fr\ 1 \^( T dadh d d\og{h a ) \ 

The tensor T appears as the dyad mv in equation (30). 

First we consider Cartesian coordinates (a: 1 , a; 2 , a: 3 ) = (x,y,z). The geometric scale factors are h x — 
h y = h z = 1 so that Sg(T) is identically zero. The cell volume is AVijk = AxiAyjAz k and the surface 
areas of the six bounding cell faces are 



A 



x i± Uk = AyjAzk , Al j±hk = A Xi Az k , A* t ._ fe± i = A Xi A yj . (A4) 

In cylindrical coordinates, (a; 1 , a; 2 , x 3 ) = (r, <f>, z), the scale factors arc 

h r = 1 , h<f, = r , h z = 1 . (A5) 

while the cell volume is 

r r i+ i r^j+i / >z fe+i / r2 ii — r2 _i\ 

AVi,j,k = 2 2 2 rdrdzd^ = l+ ~ 2 ^ 1 1 A^Az k . (A6) 

i—^ i~\ k ~i ^ ' 

The areas of cell faces are 

/r 2 - r 2 \ 

A\ ±h . k = r i±h A^Az k , Af. ±hk = A n Az k , A* . fe±| = ^ i+ * - ±i j , (A7) 
and the geometrical source vector is 

S G (mv) = ^—=-n r -\ — . (A8) 
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In spherical coordinates, (x 1 , x 2 , x 3 ) = (r,6,cj)), the scale factors are 

h r = 1 , hg = r , — r sin , (A9) 

the cell volume is 

/ r i+ i r e j+i f^fc+i /r 3 i -r 3 _i\ , . 

2 / 2 / 2 r 2 sin OdrdOdcj} = ' +2 - - 2 (cos^.i - cos^+i J Ac/) k , (A10) 

and surface areas are 

A[ ±iij>fc = r\ ±h (cos^._i - cos^+i) A0 fc , (All) 
= ( ^s^ 1 ) ™*j ±i A0 fc , (A12) 

I+ % — j A ^ • ( A1 3) 

In this case the geometrical source term is 

/ + \ „ / to iv - cot Orris v<t> \ . , f m^Vr + cot 9m (f ,vg\ ^ 
S G (mv) = I 1 n r + I \n g + I 2— I n . (A14) 

B. RECONSTRUCTION PROCEDURE 

In our parabolic reconstruction procedure we closely follow the prescription of CW84 with only minor 
modifications. The contact steepening algorithm, specific to PPM, has been used only for one-dimensional 
tests and will not be reported here (for details see Marti & Miiller 1996). We also note that, as it was pointed 
out by CW84, in order to preserve accuracy of the scheme the interpolation should be formally done using 
the conservative variables. Reconstruction based on primitive variables offers, however, some advantages 
(preserves pressure positivity and in case of relativistic hydrodynamics allows to satisfy kinematical and 
thermodynamic constraints 2 ) and is our method of choice. 

Let q be the volume-average of some quantity q, Q(£) = / <?(£') d£' its integral, where £ is a generic 
volume coordinate. Here we adopt the implied notation q <-> qi,j,k whenever i, j or k are not present (the 
same convention is adopted also for other three-dimensional quantities). The reconstruction process begins 
with interpolating a quartic polynomial through Qi± 2+ ± ,Q,±i + i >Qi+± ■ This polynomial is used to calculate 
the single-valued estimate 

dQ(0 



at the zone interfaces. The explicit result of this construction is (CW84) 



(Bl) 



= qi + <H (qi+i - qi) + M<7* - c t 5q i+1 , (B2) 



2 Specifically, by requiring that r < 2, \v\ < 1 and h > 1 one can show that E > \m\, and D/E < sf \ — m 2 /E 2 must be 
satisfied at all times. 
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with 

| 2(Ag i+lCi - A^6j) (B3) 



A& + A^+i A& + A^+i 
A^+i \ A^+i + A&+2 / A& \ A&_i + A£i 



(B4) 



Monotonicity is enforced by limiting q i+ i to lie between and (fr+i. This is achieved by using the 
limitcr (VanLeer 1997) 

_ r min(|%|,2|^ - ft.ipl^ - ^ +1 |)sgn(5gi) if (<? i+ i - - <ft_i) > , 

= S (B5) 

otherwise , 



where 



A&_! + A& + A^+i 



2A&_i+A& 2A&+1+A& _ \ 
(q i+ i - qi ) + — — — (qi - q i - 1/ 



(B6) 



A& +1 + A£ Hl ' A&_i + A& 

is the average slope. At the end of this step we initialize 

= 1i+\,R = <li+\ > ( B7 ) 

where <Zi+^ l arLC i R are the anc ^ right limiting values of q at the zone's right interface. 

Next, we constrain the interface values for zone k to lie within the extreme values found among 
all neighboring cells (Barth 1995). Construction for g i+ i L (identical procedure is followed to limit Qi-i r) 
proceeds as follows. Let <7 max and q min be, respectively, the maximum and minimum value of q~i t j,K for 
I = i - + 1, J = j - 1, j, j + 1, K = k- l,k,k+ 1. Then 

q i+hL = max (g min ,min (V^+i.l)) • (B8) 

Notice that in two and three dimensions there is a total of 9 and 27 zones involved in the limiting step. 

In the next step we limit the parabolic distribution of q to ensure that its profile remains monotonic in 
each cell. A first case when monotonicity can be violated is when q is a local maximum or minimum. In this 
case we revert to first order interpolation: 

-» q if ~q)(q- Qi+i t L) <0. (B9) 

Monotonicity can also be violated if the parabolic distribution (54) has an extremum inside the zone. 
In this case one of the interface values is adjusted so that the parabolic profile has an extremum at the other 
interface. The modified distribution also has to preserve the correct zone-average value. The final result is: 

Qi _ hR 3q- 2q i+hL if (q i+hL q t ^ R ) (q - l+ *' L - " R j > V - > , (BIO) 

q l+hL 3q-2 qi _ hR if (q i+hL - qi _ hR ) [q ' +2 ' L - < - A - I . (Bll) 

It should also be mentioned that the reconstruction algorithm may occasionally fail to respect the 
condition vf + v% + < 1. To prevent this from happening, interpolation of the velocity components reverts 
to first order whenever the total velocity exceeds the bounds provided by the neighboring cells, in a way 
similar to equation (B8). 



-27- 



B.l. Dissipation Algorithm 

Interpolation profiles are modified in presence of strong shocks to prevent unphysical oscillations. This 
is achieved by the flattening algorithm in which the interface values are modified as 

Qi+i,L =xq i +ix + ( 1 -x)q, (B12) 

Qi-^R = XQi-^,R + (^-X)Q, (B13) 

where x is a multi-dimensional flattening parameter, < x < 1. Note that for x = the accuracy of the 
method reduces to first order. The flattening parameter is computed as (Miller & Colclla 2001, 2002) 

X = mm {xi,i,j,k,X2,i,j,k,X3,i,j,K ) , (B14) 

where the minimum is taken over all /, J, K in the range i — 1 < I < i + 1, j — 1 < J < j + 1, k—1 < K < k+ 1. 
The coefficients xi, Xi an d X3 are one-dimensional. In what follows we describe the procedure only for %i; 
the remaining two values are calculated in the same way. First, we introduce a measure of the shock width 

\Pi + 2-Pi+2\ 

and shock strength 

Z X = bi+i-^-il (Bi 6) 
min(pi_i,p I+ i) 

Next we define 

Xr" = max ( 0, min f 1, /"^ ) ^ , (B17) 

so that flattening is not applied for (3 < [3 m%n . The one-dimensional flattening parameter is then restricted 
only to the regions where the flow undergoes compression and its value depends on the shock strength: 



Xi = < 



Xr-^min^ ^ ^ jj if « ll<+1 <« 1 , <) (Big) 

1 otherwise . 



For the present study we follow recommendation of Miller & Colella (2002) and adopt (3 max = 0.85, (3 min = 
0.75, Z max = 0.75, Z min = 0.25. 

Slope flattening is combined with the introduction of an explicit diffusive flux, in order to further reduce 
spurious oscillations behind strong shocks. To this purpose, the numerical fluxes (eq. [36]) in the final 
conservative update are augmented as 

F l+k - F i+ , +k v (U-U i+1 ) . (B19) 

Here 

k v = a max (-2> i+ 1 , o) , (B20) 

where a is typically set to 0.1. T> i+ 1 is a measure of the convergence of the flow at the zone interface i + In 
cartesian coordinates with a uniform zone spacing (Ax = Ay = Az), T> i+ i is a discrete undivided difference 
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approximation to the multidimensional divergence of v. In order to compute T> i+ ^ we define a measure of 
the convergence of the flow at the cell corners as 



^ {Vl,i+1,J,K - V\,i,J tK ) + ^ ( W 2,/j + l,K - V2,I,j,K) 

= k,k+l 

+ ^2 {V3,I,J,k+l - V3,I,I,k) 



7=«,»+l 
K = k,k + 1 



J=i,«+1 
J=3,j + 1 



(B2I) 



so that D i+ i is obtained by a simple average procedure 

= i (^i+JJ+J.fc+J + P i+iJ-i,fe+i + C i+i,j-i,fc-i) • 



(B22) 
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At 




t 

Fig. 1. — Solution of the Riemann problem in the £ — t plane. The discontinuity between the two states L and R initially 
located at £ = decays into a set of three waves. The emerging pattern divides the £ — t plane into the four distinct regions 
labeled by L, L* , R* , R. A contact discontinuity separates state L* from state R* , while the R* R and LL* waves can be either 
shocks or rarefactions (in this picture the LL* and R* R wave are, respectively, a rarefaction wave and a shock). 




Fig. 2. — Graphical solution in the v — p plane of v(p, R) and v(p, L); the intersection between the two curves is the solution 
to the Riemann problem in the two-shock approximation. 



-32 - 



4 



3.5 



2.5 



2 







r = 4/3 








/// Taub EOS (TM) 






/'/ / Relativistic Perfect Gas EOS (RP) 






/// Interpolated EOS (IP) 






/'/ / Ideal Gas, Constant T (ID) 






/ r = 5/3 







0.001 0.01 0.1 1 10 100 1000 

Log 8 



Fig. 3. — Plots of the function r r (0) vs. temperature logarithm (log©, © = pr) for the four different EoS discussed in the 
text. The TM EoS (solid line) marks the lower boundary of Taub's inequality (93). The relativistic perfect gas EoS (RP, 
eq. [94]) always lies above the boundary line, and therefore satisfies Taub's inequality for every 0. The interpolated EoS (IP, 
eq. [95]) approaches the correct limit only at high temperatures. The horizontal lines correspond to the ideal EoS plotted for 
T = 5/3 and T = 4/3, but only the second one satisfies Taub's inequality for < 9 < oo. 



Relativistic Shock Tube, t = 0.36 




X 



Fig. 4. — The numerical solution (given by the symbols) compared to the exact analytical solution (solid line) for the 
relativistic shock tube (PI in the text), at t=0.36 is shown. The wave structure is well described by our numerical scheme: 
both the shock and the contact wave are smeared out on 4 — 5 and 2 — 3 zones, respectively, in the density profile (asterisk 
marks). For this test problem 400 equally spaced zones were used and C a = 0.9. 
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Fig. 5. — Relativistic shock tube (P2): analytical solution (solid lines) and numerical solutions on 400 zones are shown at 
t = 0.4; from left to right v yR = 0, 0.9, 0.99, from top to bottom v yL = 0, 0.9, 0.99. 
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Fig. 6. — Initial condition for the 2-D riemann problem (on the left) and solution at t= 0.4 (right) 
simple shocks at the (2 <-» 1), (4 <-> 1) interfaces, the value of density and pressure in region 1 are pi 
pi = 2.762987 ■ 1CT 3 . 



In order to have 
= 5.477875 ■ 10" 3 , 
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Density Velocity 




Fig. 7. — Three dimensional isosurfaces of density (top left), total velocity (top right), pressure (bottom left) and specific 
internal energy (bottom right) for the relativistic spherical shock reflection test problem (RSSR) at t = 2. The inflow velocity 
for this case is Vi n = —0.9. The intensity plots below each rendered surface show two-dimensional slices of the same quantity 
in the XY plane at at z = 0. Finally, the numerical (diamonds) and analytic (solid line) solutions at y, z = are compared in 
the one dimensional plots projected on the YZ plane. The resolution adopted for this case is 101 3 and integration has been 
carried with CFL = 0.4. 
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Head Position 




Time 

Fig. 8. — Distance traveled by the jet head as a function of time. The triangles show the position of the head in our simulation 
at different times, while the solid line represents the relativistic onc-dimcnsional estimate. 



Relativistic Jgt - Lc-g{p), i = 80.00 




Fig. 9. — Density distribution, in logarithmic scale, for the axisymmctric relativistic jet at t = 80. The jet is initially in 
pressure equilibrium with the ambient medium, which is 100 times denser than the jet; the initial velocity is v z = 0.995. The 
resolution is 24 zones per beam radius. 
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Lorentz Factor, t = 3.5 s 

1 20 [ 
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Fig. 10. — Lorentz factor for the gamma ray burst problem at t = 3.5 s. The reverse shock is located at logr fs 8.9. 



y.t Lo-gtpl, 1 = -1.5 




Fig. 11. — Density distribution at t = 0.7 s (left panel) and t = 3.5 s (right panel) for the gamma ray burst problem. The 
grid is unevenly spaced with a higher resolution close to the polar axis. The low density channel along the polar axis is evident 
at both times. The result can be compared to Zhang et al. (2003), Figure 2. 
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Table 1. Expressions for the different EoS, square of sound speed and solutions to the Taub adiabat for 
the four different cases presented in the text; the coefficients a, b and c for the ID and TM EoS are given 
by equations (98) and (99). The * case has to be solved numerically. 

EoS h(Q) c 2 Taub adiabat 



RP 



h 2a 

K 3 {i/e) e i + 5he - h 2 „ 

K 2 (l/e) h l + 5hQ-h 2 -Q 2 



if 29 + v^TT h 2 = h 2 + 4 4^ 

h + 20 s 3p + Ps 



21 



4 3h h-e 11 \b\ + Vb 2 - 4ac 



Table 2. Expressions for the mass flux j, and derivative d(hr)/dp for the four different EoS presented in 
the text; here w = hr. Notice that for all but the RP case, the expressions for j 2 is numerically well 

behaved for [p] — > 0. 

EoS j 2 d(hr)/dp 



T r p w + ws- 2KhT r w 



ID 



RP 



IP 



TM 



Y r W S + (W + W S ) - l) A T rP - [p] 

pps[p] (w + w s )(h — hQ — 2hhw) 

Ps[hO}-h s Os[p} 2hhp - [p](h - hO) 

4p ws — 3u> 

3ws — w ip + ps 

\b\ + Vb 2 - Aac w(3h 2 - 3/i6 + 1) - w s {2h 2 - I - 5/i6) 



2w s {h 2 s + 2w s p + 2) 1 + [p}(l + 5hO - 2h 2 ) 
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Table 3. Discrete L-l norms of the conserved quantities for the shock tube problem (PI). 



Ax- 1 


IMi 


IMIi 


IMIi 


50 


0.1620 


0.2160 


0.1560 


100 


0.0920 


0.1420 


0.1000 


200 


0.0505 


0.0792 


0.0558 


400 


0.0298 


0.0436 


0.0319 


800 


0.0118 


0.0264 


0.0206 


1600 


0.0115 


0.0149 


0.0126 



Note. — Errors are computed accord- 
ing to ||e g ||i = Y.i Ax\qi - Q(xi)\ (as in 
Marti & Miiller (1996)), where q is one of 
D, m, E, and Q(xi) is the exact solution 



Table 4. Relative global errors (computed as in Aloy et al. 1999) for the Relativistic Spherical Shock 
Reflection (RSSR) test problem for different inflow velocities at t = 2. Computations have been performed 

on a 81 3 grid. 



v 7 IMIi IMIi IMIi 



10 


-1 


2.29 


0.18 


0.01 


0.29 


10 


-3 


22.37 


0.20 


0.03 


0.17 


10 


-5 


223.6 


0.26 


0.04 


0.20 


10 


-7 


2236.1 


0.21 


0.04 


0.24 



Note. - The inflow velocity is given by 
Vin = v — \. The corresponding Lorcntz factor is 
given in the second column. The first two cases 
have been run with CFL = 0.4, while CFL = 0.1 
has been used for the last two cases. 



